
# ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
# DESCRIPTION
# ______________________________________________________________________________

# Author: Changwook Ju, Yale University (changwook.ju@yale.edu)

# This code creates Figure 1.

# Preamble ---------------------------------------------------------------------

rm(list = ls())

# Packages ---------------------------------------------------------------------

library(tidyverse)
library(ggh4x)
library(ggtext)
library(haven)
library(patchwork)

#### function to make center axis title with double plots

add_global_label <-
  function(pwobj,
           Xlab = NULL,
           Ylab = NULL,
           Xgap = 0.03,
           Ygap = 0.03,
           ...) {
    ylabgrob <- patchwork::plot_spacer()
    if (!is.null(Ylab)) {
      ylabgrob <- ggplot() +
        geom_text(aes(x = .5, y = .5), label = Ylab, angle = 90, ...) +
        theme_void()
    }
    if (!is.null(Xlab)) {
      xlabgrob <- ggplot() +
        geom_text(aes(x = .5, y = .5), label = Xlab, ...) +
        theme_void()
    }
    if (!is.null(Ylab) & is.null(Xlab)) {
      return((ylabgrob + patchworkGrob(pwobj)) +
               patchwork::plot_layout(widths = 100 * c(Ygap, 1 - Ygap)))
    }
    if (is.null(Ylab) & !is.null(Xlab)) {
      return((ylabgrob + pwobj) +
               (xlabgrob) +
               patchwork::plot_layout(
                 heights = 100 * c(1 - Xgap, Xgap),
                 widths = c(0, 100),
                 design = "
                                   AB
                                   CC
                                   "
               )
      )
    }
    if (!is.null(Ylab) & !is.null(Xlab)) {
      return((ylabgrob + pwobj) +
               (xlabgrob) +
               patchwork::plot_layout(
                 heights = 100 * c(1 - Xgap, Xgap),
                 widths = 100 * c(Ygap, 1 - Ygap),
                 design = "
                                   AB
                                   CC
                                   "
               )
      )
    }
    return(pwobj)
  }

# Data Import ------------------------------------------------------------------

data <- read_dta("data/Faulkner_Welsh_2022.dta")

# Plot -------------------------------------------------------------------------

## US DoS ----------------------------------------------------------------------

data$facet <- "US DoS"

label_prev <- c("0", "1", "2", "3")

## violin

violin_statedept_sv <- data %>%
  filter(!is.na(statedept_sv)) %>%
  ggplot(aes(x = factor(statedept_sv))) +
  geom_violin(aes(y = year, fill = factor(statedept_sv)),
              scale = "width",
              linewidth = 0.5) +
  geom_jitter(
    aes(y = year),
    alpha = 0.2,
    color = "black",
    fill = NA,
    shape = 1,
    width = 0.15
  ) +
  geom_boxplot(
    aes(y = year),
    color = "#f2f2f2",
    fill = NA,
    alpha = 0.8,
    width = 0.3,
    outlier.shape = NA
  ) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  facet_wrap(. ~ facet) +
  scale_x_discrete(labels = label_prev) +
  labs(y = "Year") +
  theme_bw() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    strip.background = element_rect(
      linewidth = 0.5,
      fill = "white",
      color = "black"
    ),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_text(size = 16),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.line = element_line(),
    axis.text.x = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 16),
    text = element_text(size = 16),
    axis.title.y = element_text(margin = margin(
      t = 0,
      r = 10,
      b = 0,
      l = 0
    ))
  ) +
  coord_cartesian(clip = "off")

## bar

bar_statedept_sv <- data %>%
  filter(!is.na(statedept_sv)) %>%
  ggplot(aes(x = factor(statedept_sv))) +
  geom_bar(
    aes(fill = factor(statedept_sv)),
    linewidth = 0.5,
    width = 0.8,
    color = "black",
    show.legend = FALSE
  ) +
  geom_richtext(
    stat = "count",
    size = 5,
    vjust = 0.1,
    lineheight = 0.8,
    fill = NA,
    label.color = NA,
    aes(
      label = paste0(
        "<span style='font-size:13pt'>",
        scales::percent(after_stat(count) / sum(after_stat(count)),
                        accuracy = 0.1),
        "</span>",
        "<br>",
        after_stat(count)
      )
    )
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 330)) +
  scale_x_discrete(labels = label_prev) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  labs(x = "Prevalence of Rebel CRSV",
       y = "Count") +
  theme_bw() +
  theme(
    strip.background = element_rect(color = NA, fill = "white"),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_blank(),
    axis.line = element_line(),
    text = element_text(size = 16)
  ) +
  coord_cartesian(clip = "off")

## AI --------------------------------------------------------------------------

data$facet <- "AI"

label_prev <- c("0", "1", "2", "3")

## violin

violin_ai_sv <- data %>%
  filter(!is.na(ai_sv)) %>%
  ggplot(aes(x = factor(ai_sv))) +
  geom_violin(aes(y = year, fill = factor(ai_sv)),
              scale = "width",
              linewidth = 0.5) +
  geom_jitter(
    aes(y = year),
    alpha = 0.2,
    color = "black",
    fill = NA,
    shape = 1,
    width = 0.15
  ) +
  geom_boxplot(
    aes(y = year),
    color = "#f2f2f2",
    fill = NA,
    alpha = 0.8,
    width = 0.3,
    outlier.shape = NA
  ) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  facet_wrap(. ~ facet) +
  scale_x_discrete(labels = label_prev) +
  labs(y = "Year") +
  theme_bw() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    strip.background = element_rect(
      linewidth = 0.5,
      fill = "white",
      color = "black"
    ),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_text(size = 16),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.line = element_line(),
    axis.text.x = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 16),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

## bar

bar_ai_sv <- data %>%
  filter(!is.na(ai_sv)) %>%
  ggplot(aes(x = factor(ai_sv))) +
  geom_bar(
    aes(fill = factor(ai_sv)),
    linewidth = 0.5,
    width = 0.8,
    color = "black",
    show.legend = FALSE
  ) +
  geom_richtext(
    stat = "count",
    size = 5,
    vjust = 0.1,
    lineheight = 0.8,
    fill = NA,
    label.color = NA,
    aes(
      label = paste0(
        "<span style='font-size:13pt'>",
        scales::percent(after_stat(count) / sum(after_stat(count)),
                        accuracy = 0.1),
        "</span>",
        "<br>",
        after_stat(count)
      )
    )
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 330)) +
  scale_x_discrete(labels = label_prev) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  labs(x = "Prevalence of Rebel CRSV",
       y = "Count") +
  theme_bw() +
  theme(
    strip.background = element_rect(color = NA, fill = "white"),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_blank(),
    axis.line = element_line(),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

## HRW -------------------------------------------------------------------------

data$facet <- "HRW"

label_prev <- c("0", "1", "2", "3")

## violin

violin_hrw_sv <- data %>%
  filter(!is.na(hrw_sv)) %>%
  ggplot(aes(x = factor(hrw_sv))) +
  geom_violin(aes(y = year, fill = factor(hrw_sv)),
              scale = "width",
              linewidth = 0.5) +
  geom_jitter(
    aes(y = year),
    alpha = 0.2,
    color = "black",
    fill = NA,
    shape = 1,
    width = 0.15
  ) +
  geom_boxplot(
    aes(y = year),
    color = "#f2f2f2",
    fill = NA,
    alpha = 0.8,
    width = 0.3,
    outlier.shape = NA
  ) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  facet_wrap(. ~ facet) +
  scale_x_discrete(labels = label_prev) +
  labs(y = "Year") +
  theme_bw() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    strip.background = element_rect(
      linewidth = 0.5,
      fill = "white",
      color = "black"
    ),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_text(size = 16),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.line = element_line(),
    axis.text.x = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 16),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

## bar

bar_hrw_sv <- data %>%
  filter(!is.na(hrw_sv)) %>%
  ggplot(aes(x = factor(hrw_sv))) +
  geom_bar(
    aes(fill = factor(hrw_sv)),
    linewidth = 0.5,
    width = 0.8,
    color = "black",
    show.legend = FALSE
  ) +
  geom_richtext(
    stat = "count",
    size = 5,
    vjust = 0.1,
    lineheight = 0.8,
    fill = NA,
    label.color = NA,
    aes(
      label = paste0(
        "<span style='font-size:13pt'>",
        scales::percent(after_stat(count) / sum(after_stat(count)),
                        accuracy = 0.1),
        "</span>",
        "<br>",
        after_stat(count)
      )
    )
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 330)) +
  scale_x_discrete(labels = label_prev) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  labs(x = "Prevalence of Rebel CRSV",
       y = "Count") +
  theme_bw() +
  theme(
    strip.background = element_rect(color = NA, fill = "white"),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_blank(),
    axis.line = element_line(),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

## Max -------------------------------------------------------------------------

data$max_sv <- pmax(data$statedept_sv,
                    data$ai_sv,
                    data$hrw_sv) # create `max_sv`

data$facet <- "Max"

label_prev <- c("0", "1", "2", "3")

## violin

violin_max_sv <- data %>%
  filter(!is.na(max_sv)) %>%
  ggplot(aes(x = factor(max_sv))) +
  geom_violin(aes(y = year, fill = factor(max_sv)),
              scale = "width",
              linewidth = 0.5) +
  geom_jitter(
    aes(y = year),
    alpha = 0.2,
    color = "black",
    fill = NA,
    shape = 1,
    width = 0.15
  ) +
  geom_boxplot(
    aes(y = year),
    color = "#f2f2f2",
    fill = NA,
    alpha = 0.8,
    width = 0.3,
    outlier.shape = NA
  ) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  facet_wrap(. ~ facet) +
  scale_x_discrete(labels = label_prev) +
  labs(y = "Year") +
  theme_bw() +
  theme(
    panel.border = element_rect(fill = NA, color = "black"),
    strip.background = element_rect(
      linewidth = 0.5,
      fill = "white",
      color = "black"
    ),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_text(size = 16),
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.line = element_line(),
    axis.text.x = element_blank(),
    plot.title = element_text(hjust = 0.5, size = 16),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

## bar

bar_max_sv <- data %>%
  filter(!is.na(max_sv)) %>%
  ggplot(aes(x = factor(max_sv))) +
  geom_bar(
    aes(fill = factor(max_sv)),
    linewidth = 0.5,
    width = 0.8,
    color = "black",
    show.legend = FALSE
  ) +
  geom_richtext(
    stat = "count",
    size = 5,
    vjust = 0.1,
    lineheight = 0.8,
    fill = NA,
    label.color = NA,
    aes(
      label = paste0(
        "<span style='font-size:13pt'>",
        scales::percent(after_stat(count) / sum(after_stat(count)),
                        accuracy = 0.1),
        "</span>",
        "<br>",
        after_stat(count)
      )
    )
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 330)) +
  scale_x_discrete(labels = label_prev) +
  scale_fill_manual(values = c(
    "0" = "#B8D9A9FF",
    "1" = "#8DBC80FF",
    "2" = "#5D9D52FF",
    "3" = "#287A22FF"
  )) +
  labs(x = "Prevalence of Rebel CRSV",
       y = "Count") +
  theme_bw() +
  theme(
    strip.background = element_rect(color = NA, fill = "white"),
    axis.ticks = element_line(color = "black"),
    axis.text = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    strip.text = element_blank(),
    axis.line = element_line(),
    text = element_text(size = 16),
    axis.title.y = element_blank()
  ) +
  coord_cartesian(clip = "off")

# Patchwork --------------------------------------------------------------------

(violin_statedept_sv / bar_statedept_sv) |
  (violin_ai_sv / bar_ai_sv) |
  (violin_hrw_sv / bar_hrw_sv) |
  (violin_max_sv / bar_max_sv) &
  theme(plot.margin = margin(5, 5, 1, 4, "pt"))

ggsave(
  "figures/pirate plots/pirate_plot_composite.pdf",
  device = cairo_pdf,
  width = 14,
  height = 7
)
